BMC Bioinformatics

Motivation: Real-time qPCR data are commonly analyzed to obtain a Cq values, which are used to compare the relative abundance of a targeted template. However, improper background subtraction distorts Cq and changing the assigned threshold can alter experimental conclusions. In addition, a Cq values provide no information on the quality of the qPCR for a given sample. Alternatively, global fitting of qPCR data provides agnostic quantification without the use of Cq or thresholds, and the fitting terms provide a useful metric of reaction performance. Current implementations of global fitting are limited to custom in-house programs. Results: qPyCR is comprehensive analysis tool written in Python that performs background signal subtractions, computes a threshold, and reports Cq values in a traditional manner. The program also performs global fitting and reports template abundances along with reaction quality values. The quality values can be used to assign performance standards, which allow researchers to distinguish aberrant reactions from true changes in template abundance. Availability and Implementation: Code and example data are freely available on GitHub (https://github.com/sdmoore-labs/qPyCR) and on Zenodo (https://doi.org/10.5281/zenodo.21115342).

BackgroundSecondary haplotigs retained in a primary assembly inflate genome size and duplicate gene content, and their removal remains a necessary curation step for heterozygous diploid genomes. HapSolo (Solares et al. 2021) treats this step as an optimization problem: it searches three pairwise alignment thresholds (percent identity, query coverage, and the ratio of aligned query to aligned reference length) to minimize a cost defined over single-copy, duplicated, and missing conserved orthologs. The original implementation depended on a separate BUSCO v3 run for every contig, on Python 2.7 and pandas in its inner loop, on a single randomized forward-walk optimizer, and on OrthoDB v9 lineage datasets that are no longer distributed. ResultsHere we report HapSolo2, a reimplementation that replaces the per-contig BUSCO step with miniprot-based ortholog classification against OrthoDB v9, v10, or v12.2 lineage datasets, adds steepest descent and simulated annealing to the original random forward walk, and evaluates all walkers in batched CuPy array operations on a single GPU. On an avocado assembly of 5,122 contigs scored against 1,440 orthologs, 100,000 iterations of random walk and steepest descent ran 5.6- and 7.6-fold faster on one NVIDIA RTX A6000 than on 10 CPU processes. Across 45 runs on four diploid assemblies (Anopheles funestus, Persea americana, Vitis vinifera, and Amblyraja radi-ata; 212 Mb to 3.2 Gb) and three OrthoDB releases, the three optimizers reached best costs within 1.8% of one another. Independent assessment with BUSCO v6.1 showed that HapSolo2 reduced duplicated orthologs from 52.0% to 5.6% in P. americana, from 17.9% to 10.4% in V. vinifera, and from 7.8% to 2.8% in A. funestus, with completeness decreasing by at most 1.3 percentage points. On the avocado assembly, HapSolo2 reduced the cost reached by the original HapSolo from 0.2200 to 0.1036 and raised contig N50 from 3,366 kb to 4,608 kb. ConclusionsHapSolo2 makes ortholog-guided haplotig removal practical in minutes on a single workstation, independent of a BUSCO installation, and reproducible across OrthoDB releases. It is distributed as an installable Python package with a unified command-line pipeline, an assembly statistics and plotting tool, an Apptainer container, and archived copies of the legacy ortholog datasets.

We have developed a method to integrate transcriptomic coexpression networks across diverse experimental conditions within a single species. Our framework expands previous pannetwork approaches into a hypergraph based pannetwork. It first identifies coexpressed gene clusters within each individual dataset and extracts frequently co-expressed gene sets across multiple datasets using a frequent itemset mining algorithm. This process yields a hypergraph where each hyperedge is assigned a frequency (universality, U). We then extracted a subnetwork comprising high U hyperedges as the core network and formed modules within it. We applied our method to 106 Escherichia coli transcriptome datasets from the GEO database. The modularity of the core network peaked at a universality cutoff of 15, which was subsequently used to define it. Approximately 70% of the resulting core modules correspond to operons, and conversely, approximately 70% of all operons are covered by these core modules. We visualized the core network via an inter-modular network and analyzed core module dataset relationships using a modularity profile matrix. Based on these analyses, we successfully visualized the dynamic reorganization of the coexpression network in response to environmental changes in bacteria.

In recent years, the issue of reproducibility in scientific studies has been at the forefront of scientific discussion. Many of these issues arise due to poor employment of statistical methods, making it crucial to address these knowledge gaps for frequently used methods. Elastic net regression (ENET) is one such algorithm used frequently in biology and public health due to its ability to select important predictors (also known as "features") from large data sets. However, ENET struggles with consistent feature selection and hyperparameter choice in these applications due to low observation counts and noisy data. This limitation hinders ENET's effectiveness in biological research due to low model reproducibility and predictor interpretability. We propose the Instability Quotient Criterion (IQC) as a metric to quantify the consistency of selected features for feature selection models such as ENET, with a specific emphasis on biological applications. IQC leverages ensemble modeling and Shannon's entropy in order to determine the stability of selected features derived from ENET models. We demonstrate IQC's robustness via simulation studies, showing that IQC is capable of detecting significant decreases in model stability with lower observation counts, higher dimensionality, and higher proportions of noisy signal. We additionally utilize simulation to provide a naive threshold for low versus high stability. Finally, we applied IQC to a case study to illustrate how our method can be integrated into a larger analysis pipeline when looking at microarray gene expression data. IQC indicated stability issues for differentially-expressed genes when classifying acute leukemia subtypes in spite of high classification success. Lastly, we discuss how model instability impedes interpretation of biologically significant genes for discerning cancer subtypes. Overall, we present IQC as a robust evaluation metric that can easily be incorporated into high-dimensional biological analyses.

Genetic association analyses of exome sequencing data pose distinct methodological challenges in small cohorts, which are common in rare disease research. Data filtering prior to association testing requires particular care, especially when analyzing datasets generated under different study designs. Most of current genetic association tools need to be integrated into complex data processing pipelines, whose implementation requires technical expertise and is prone to error. To alleviate mentioned issues we built a comprehensive, high performance Nextflow workflow for genetic association testing, which includes all commonly performed data filtering and preparation steps, controlled by a rich set of modifiable parameters. We compared our software against existing tools incorporating quality control and aggregated association testing steps and found our solution to offer the most complete suite of data preparation, association testing and reporting modules, making small-cohort association analyses accessible to teams with limited technical expertise.

Genomic variant calling is a central computational task in disease research and diagnostics, where somatic mutations such as JAK2 V617F serve as diagnostic markers in myeloid neoplasms. VarDict, originally written in Perl, and its reimplementation VarDictJava are widely used variant callers in the bioinformatics community. The Perl implementation requires run times of up to several hours; VarDictJava reduces this to a few minutes but can consume large amounts of main memory, which limits the number of jobs that can be run in parallel on a given machine. We present vardictcpp, a functionally identical reimplementation of VarDictJava in C++17. vardictcpp is approximately four times faster in single-core execution and 9-11 times faster using eight cores, while requiring up to 17 times less main memory. Variant calls are identical to those of VarDictJava on the evaluated whole-exome samples. The reduced memory footprint allows more concurrent jobs on existing hardware and thereby increases sample throughput, which is particularly relevant for clinical sequencing laboratories with fixed compute resources. vardictcpp is implemented in C++17 and freely available at https://github.com/MHH-Bioinformatics-Hematology/vardictcpp

Several tools are available to detect differential gene expression. The predominant approach of different statistical methods is to consider changes in gene expression independent from their genomic context. Previously, DELocal was introduced that takes into consideration gene expression dynamics in the genomic neighbourhood. Whereas DELocal was shown to identify developmentally regulated genes in developing mouse tooth, its performance in other biological conditions remains undetermined. To benchmark DElocal together with other methods, we develop another tool 'broadSeq'. It provides easy to use interface for methods like DELocal , deseq2, limma, edger, EBSeq and many more. Here we report benchmarking results on different data sets that comprise wild type developing mouse organs from different embryonic time points, and also human cancer data. The examined datasets are from the heart, kidney, liver, brain, forebrain, and hindbrain. We found that most methods produce relatively comparable lists of differentially expressed genes. Nevertheless, DELocal identifies partially distinct sets of genes compared to the other methods. Evaluating the identified genes using gene and disease ontology terms show that the uniquely detected genes by DELocal are biologically relevant and also highly expressed compared to the genes identified by the other methods. Because consistently highly expressed genes show relatively subtle fold changes, the neighbourhood concept of DELocal appears to correct for the difficulty in detecting these differential expression of highly expressed genes. Finally, we provide the tools as combined bioconductor packages, allowing researchers to explore and decide the best methods for their data and questions.

The advent of high-throughput sequencing technologies has made it essential to employ advanced tools for data integration and interpretation. In this work, we propose to build upon two existing efforts by expanding their scope and applicability to multi-omic data integration through the development of a pipeline that generates a topologically informed ranking of genes to assess gene relevance to disease. Specifically, we exploit topological data analysis, capitalizing on: its graph-based framework, to refine complex interactions characteristic of biological systems; its intrinsic explainability, to enable interpretation of model outputs in terms of biologically meaningful features; its robustness to small input perturbations, to ensure stable results in the presence of noisy data, as omics are. In particular, we employ Harmonic Persistent Homology (HPH), and propose two relevant advances. First, we adopt an integration and clustering tool (i.e. iNETgrate) as an informed means to reduce data size and integrate multiple omic layers by aggregating methylation beta values from loci to the gene level, to be combined with gene expression values. This enables the application of HPH, otherwise limited by the computational burden. Second, as a non-trivial consequence, we apply HPH to multi-omic molecular data, and not to patients as done so far, enlarging HPH scope. Finally, we follow up on recent efforts toward performance standardization, and validate our results against biomarkers from a well-known TCGA breast cancer benchmark dataset.

Multiple sequence alignments (MSAs) are often treated as a solved problem. However, ambiguously or incorrectly aligned regions are the source of much error in phylogenetics and comparative genomics. CRAIC is a cross-platform desktop workbench that treats an alignment as a set of hypotheses to be stress-tested. CRAIC uses the posterior match probabilities of a three-state pair hidden Markov model (pair-HMM) to derive complementary views of alignment confidence. The five views are (i) a reliability overlay with live column and residue masking; (ii) a multi-aligner disagreement map, (iii) a local realignment sandbox, (iv) an interactive posterior explorer, and, where the true alignment is known, (v) a per-residue view of correctly predicted homologies. When present, MAFFT, MUSCLE, ProbCons, Clustal Omega, ClustalW and PRANK are integrated in the tool. In benchmarks the reliability score can locate misaligned columns, but while CRAIC can implement automatic filtering, it is primarily built for inspection and curation.

Bioinformatics software plays a central role in biological research, yet many academic projects lack robust testing and continuous integration (CI) practices. This absence can compromise reproducibility, accessibility, and long-term maintainability. In this study, we conduct an empirical analysis of GitHub repositories spanning bioinformatics and other scientific software domains to examine the relationship between CI, software testing, and software quality. Repository discovery was performed using domain-specific search terms via GitHub's Application Programming Interface (API). This was followed by a layered data extraction pipeline, which included a Large Language Model (LLM) validation step, that connected repositories to published works, institutions, authors, funding, and citation counts. By integrating repository activity, community engagement, and backlog health into a composite Repository Quality Index (RQI), we find that the presence of CI and software tests is significantly associated with higher repository quality. Additionally, repositories with higher academic citation counts tend to cluster in higher RQI tiers, suggesting that software quality, which is driven by proper testing and integration, is associated with academic impact and funding status. These results demonstrate that adopting industry-standard development practices, such as CI and software testing, is a tractable path toward more sustainable and impactful academic software and funding opportunities.

Motivation: High-throughput sequencing technologies have driven rapid growth of biological sequence databases. Public repositories must therefore rely on automated computational heuristics to screen submitted sequences for errors and low quality. For example, SILVA SSU Ref, which exploits the conserved nature of 16S and 18S rRNA sequences, applies strict algorithmic quality controls yet still accepts sequences with up to 30% of their nucleotides deviating from any previously accepted sequence. This permissiveness creates opportunities for the admission of modified sequences, such as biologically plausible sequences generated by DNA foundation models. The vulnerability of public sequence databases to becoming polluted or poisoned with modified sequences necessitates the development of methods to detect such sequences. Results: We present the first investigation, to our knowledge, of the detectability of modified sequences. We consider simple computationally-modified 16S rRNA sequences that pass the quality control inclusion criteria of the SILVA database, which we generate via random substitutions. We present classifiers that can distinguish such modified sequences from natural 16S rRNA using conserved motifs. Our best classifier achieves over 90% sensitivity and specificity on our testing set when using a 5% artificial mutation rate. One feature used in our classifiers, gapped k-mers constructed from universally conserved nucleotides, was conserved across all three domains of life despite relying on exact matches to patterns found in E. coli, advancing our understanding of conserved grammatical structure in small subunit rRNA sequences. Availability and Implementation: Our source code is available at https://github.com/rainhaworth/16S-Mutation-Classifiers.

Phenotype prediction from public gene-expression data is straightforward to implement but reliable prediction is harder. Leakage during feature selection can inflate performance estimates. Class imbalance can hide poor performance in minority classes. Confounds like tissue type, repeated sampling, and demographic differences can also create patterns that resemble phenotype-associated biology. Here, I present LOL (Latent Omics Learning), a Python package for phenotype prediction from bulk GEO RNA-seq and microarray data. LOL uses leakage-safe preprocessing, imbalance-aware evaluation, and group-aware validation where repeated samples are present. Each analysis produces two audit diagnostics. The leakage inflation index ({Delta}LII) measures the change in performance between naive and leakage-safe evaluation. PVCA-lite examines associations between the latent representation and available metadata covariates, including phenotype. I evaluated the workflow across six independent GEO cohorts. The datasets covered microarray and RNA-seq platforms, binary and multi-class phenotypes, and sample sizes ranging from 18 to 566. The analyses produced a varied set of outcomes. Across diverse GEO datasets, LOL can reveal phenotype-associated expression patterns while identifying factors distort their interpretation. These include information leakage, class imbalance, repeated sampling, tissue differences, and demographic or other metadata-associated confounding. The resulting diagnostics help distinguish robust predictions from results that require further scrutiny. LOL does not introduce a new representation-learning algorithm. Its contribution is a reproducible workflow that places validation and diagnostic checks alongside prediction. The resulting analysis can expose weaknesses in an apparently convincing result and make those limitations visible. LOL is freely available under the MIT license from https://github.com/ngangao/lol-omics.

Abstract Accurate taxonomic classification in metagenomics is often compromised by low-complexity sequences, which lead to chance matches that in turn cause sequences to be misclassified. Here we present Kmask, an entropy-based masking tool implemented for use either standalone or as part of Kraken [1,2] database construction, which replaces low-entropy regions with Ns. Using a sliding window size aligned with Kraken's default k-mer length and parameters optimized across 12 control bacterial genomes spanning a broad range of GC content, Kmask efficiently removes low-complexity sequences while retaining high-complexity regions. To benchmark performance, we applied Kmask to a newly constructed database, Microbial2025, that contains over 71,000 bacterial, archaeal, viral, and fungal genomes, and we then classified human reads against both masked and unmasked versions of the database using KrakenUniq [2]. We found that Kmask substantially reduced misclassifications, driving down the false positive rate to 5.78% from 7.52%. Notably, Kmask performed comparably to an SDUST-masked [3] database, achieving a similar false positive rate (5.78% vs. 5.17%) while masking out fewer bases (1.33% vs. 1.85%). We also tested Kmask on a database of human cancer sequences, where we found that it eliminated many false positives caused by low-complexity matches between bacterial genomes and human DNA. These results demonstrate that Kmask is an effective method for masking low-complexity sequences in large microbial databases, thus improving the accuracy of metagenomic classification.

Background: Pangenome variations graphs (PVG) are promising structures to study diversity and reduce reference bias in genomics. However, it is not currently possible include complete annotation of the PVG. In a previous work, we introduced GrAnnoT, a tool for transferring the annotation from a reference genome to a PVG, and from a PVG to its embedded genomes. However, GrAnnoT was limited to transferring one annotation to the graph, and thus could not include all genes. Results: We implemented new options to GrAnnoT, that can now handle annotations of multiple genomes, and integrate them in the graph. We also developed GrAPEx, a pipeline that identifies, extracts and annotates non-annotated regions in the graph. We tested those two tools to attest their performances. Our benchmarks revealed that both GrAnnoT and GrAPEx are reliable. GrAnnoT is now able to integrate correctly several annotations, increasing the amount of biological information available in the graph. Additionally, it offers useful options to study gene distribution in the PVG. GrAPEx allows rapid annotation of the whole PVG, with no reduction in the annotation quality. Conclusion: Overall, the two tools presented here allow to annotate thoroughly the whole PVG, and integrate biological annotation to the pangenome. It also reduces the previous reference bias, where only one annotation (often the reference one) could be included in the graph. This constitutes a precious addition to the pangenomic toolkit.

Research in bioinformatics demands that many workflows and operations be consistently repeated across different tools. These operations include managing a taxdump, dealing with accession mappings, and performing analysis on the taxonomic tree. We introduce taxutils to address and alleviate the strain of these common workflows. To exhibit the extent of the package's capabilities, we introduce a strategy for assessing quality of an assembly via its potential impact on and perturbation of its mapping to a lowest-common-ancestor index. We illustrate the steps required for such an analysis, and show how the tool makes data management easy and otherwise-complex workflows simple.

In network medicine, the protein-protein interaction network, or interactome, is an essential resource for identifying candidate proteins underlying diseases and other phenotypic traits. Indeed, it can be leveraged via the guilt-by-association paradigm, which asserts that interacting proteins are likely to participate in the same molecular processes. Network propagation that combines the topology of the interactome with prior knowledge about disease genes is a promising strategy to identify new candidate genes contributing to diseases. However, existing network propagation algorithms are often biased toward highly connected proteins, called "hubs". We present BFWalk, a novel network propagation algorithm designed to avoid inflated scores for hubs. We tested BFWalk across four human phenotypes: multiple morphological abnormalities of the sperm flagella, dyschromatopsia, hypertrophic cardiomyopathy, and chronic kidney disease. For each phenotype, we assessed whether BFWalk could recover known phenotype-associated genes and whether new candidate genes were enriched in relevant tissues. Depending on the phenotype, BFWalk outperformed or matched the state-of-the-art network propagation methods. Furthermore, the results confirmed that BFWalk is free of bias toward high-degree nodes. Therefore, BFWalk is a powerful method for identifying new genes across diverse phenotypes and constitutes a strong alternative to existing network propagation algorithms. Availability: BFWalk is implemented in Python and C. The code is available under the GNU GPL (v3.0) on: https://github.com/jedrzejkubica/BFWalk

Dysfunctions in transcription cause severe human diseases. In eukaryotes, transcription initiation by RNA polymerase (Pol) II requires the formation of the pre-initiation complex (PIC) composed of general transcription factors (GTFs), Pol II, and an essential multiprotein coactivator, Mediator. In vitro experiments with GTFs and Pol II have depicted a linear sequence of PIC assembly. However, PIC assembly paths in vivo at the genomic scale and Mediator's role in these paths remain poorly understood. One strategy characterizes ChIP-seq occupancy alterations in Mediator mutants, but analysis of these data still lacks a formalized interpretation framework. In this study, we address this need by developing a quantitative framework that scores candidate assembly paths according to their ability to explain these alterations. This framework combines a model of PIC assembly based on a system of differential equations with parameter optimization based on a penalized criterion. We applied this new framework to datasets from two yeast Mediator mutants, in the Med17 and Med10 subunits. Our results show that the best-supported PIC assembly paths differ from those inferred in vitro. In particular, we show that it is important to distinguish the core TFIIH module from the TFIIK kinase module, and that similar paths starting with Mediator or TBP receive different levels of support. Furthermore, although multiple paths fit the data equally well, they differ markedly for specific gene groups. Taken together, the results show that the path can depend on the gene, and are also compatible with more than one path being operational in vivo for a single gene. This work improves our understanding of PIC assembly paths and transcription initiation, and opens perspectives for future modeling and experimental studies.

Microbial sequencing data are subject to sample-specific and species-specific distortions that complicate inference about true microbial abundances. Recent work established assumptions under which relative fold-differences in microbial abundances can be estimated from sequencing data, accompanied by an estimation and testing method called radEmu. While radEmu has excellent error rates, its computational cost can be prohibitive for large-scale analyses. Here we introduce fastEmu, a scalable approximation to radEmu's testing methodology. fastEmu uses a small, data-driven subset of taxa to approximate the computationally expensive test statistic. In a metagenomic analysis, we show a reduction in per-taxon runtime from 4-8 hours to 0.5-45 seconds, and simulations illustrate no adverse impact on estimation error or Type I error rate control. fastEmu enables testing differential abundance for thousands of microbial taxa while retaining clear parameter interpretation and error rate control.

This study introduces a novel model that integrates the capacity of Potts models to capture structural dependencies with the computational efficiency afforded by Felsenstein's algorithm in phylogenetics. While the present work focuses on protein models, the PhyloPotts framework is equally applicable to alignments of discrete states, including RNA, DNA, and 3DI characters. Our findings indicate that, although phylogenetic estimates do not consistently exhibit reduced bias relative to standard models, the approach substantially enhances the accuracy of ancestral sequence reconstruction. The method is freely available via an open source implementation in the PhyloPotts package for BEAST 2 from https://github.com/rbouckaert/potts.

Recent biological research has shown significant interest and progress towards adoption of agentic AI systems. Domain specific agentic platforms with access to large collections of biological tools with code generation and/or querying capabilities, such as Biomni (illustrative case in point), are representative of the capability to orchestrate automated research workflows. Such platforms enable collaborative AI workspace for drug discovery, genomics, bioinformatics, literature research, and experimental planning for biologists and researchers. When such platforms are used for the aforementioned workflows, individually validating every agentic decision and the research process by which it was reached becomes challenging due to the scale of decisions involved and oftentimes, insufficient disclosure of assumptions and methodological choices, limited traceability of evidence, and interpretive judgments embedded throughout the workflow, giving rise to verification debt. Verification debt in agentic environments could result in unexpected failures and consequences in downstream experiments. This study focuses on empirical qualitative characterization of verification debt in agentic (prompt to outcomes) workflows. Biomni was chosen the testing AI agentic platform in this study mainly due to the descriptive traces it provides to track every agentic step. Verification debt is examined through two complementary angles: the individual agentic decision steps that carry claims and evidence, and the underlying research process by which the agent arrives at and reports them. This study attempts to characterize verification debt with 10 test cases (covering proteomics, genomics and drug discovery) with varying degree of biological nuance across three identical and independent runs. Verification debt is represented through 9 decision-step metrics and 13 research-process metrics, and also propose the VERIFY framework for evaluating scientific verifiability and validity in agentic workflows.